clear all;
clc;
tstart=1; %start time
tmax=250; %max time of simulation
time=tstart:tmax; %time vector

n_iterations=1;
%nparticles = 1000;
nparticles = 1;
npos=2*length(time)+3; %size of pos arrays = 2t+3
pos=1:npos;
nmid_rel=tmax+2; %relative mid position in the below a arrays which corresponds to zero position in the absolute sense
nstart_rel=0;%relative start position in the below a arrays which corresponds to least -ve position in the absolute sense
nend_rel=0;%relative end position in the below a arrays which corresponds to max +ve position in the absolute sense

a_n_u_prev=zeros(npos); %a's for spin up at each N at previous time
a_n_d_prev=zeros(npos);%a's for spin down at each N at previous time
a_n_u_cur=zeros(npos);%a's for spin up at each N at current time
a_n_d_cur=zeros(npos);%a's for spin down at each N at current time

%sigma=zeros(1,npos); %the standard deviation of the spatial distribution
sigma=zeros(1,length(time));
sigma_up=zeros(1,length(time));
sigma_down=zeros(1,length(time));

avg_pos=zeros(1,length(time));
avg_pos_up=zeros(1,length(time));
avg_pos_down=zeros(1,length(time));

sigma_count = 1;
% T=tmax/2:100:tmax;
T=time;
pos_counter=0;
a=0;
epsilon = 0.05;
p_n_tmax = zeros(nparticles,npos);

p_n_u_tmax = zeros(nparticles,npos);
p_n_d_tmax = zeros(nparticles,npos);

p_n = zeros(nparticles,npos);
p_n_u= zeros(nparticles,npos);
p_n_d= zeros(nparticles,npos);

sig=zeros(1,length(time));
sig_up=zeros(1,length(time));
sig_down=zeros(1,length(time));

%for all particles
for p=1:nparticles
    for iter=1:n_iterations
        %form the initial a vectors for time = 0. consider all the particles to be
        %with spin up
    %     a = 1/2;
        a = 1/(sqrt(2));
        a_n_u_prev(nmid_rel)=a;
        a_n_d_prev(nmid_rel)=(a)*i;

        %for all times
        for t=time
            display(t);
            pos_counter = 0;
            nstart_rel = nmid_rel - t;
            nend_rel = nmid_rel + t;
            if(t ~= 1)
                a=(1/(sqrt(2))) + epsilon*rand;
            end;

            %for all possible positions for this particular time, fill up the a
            %arrays for the current time

            for pos_counter=nstart_rel:nend_rel

                a_n_u_cur(pos_counter)=  ( sqrt(1-a*a) * a_n_d_prev(pos_counter-1) )+ ( a * a_n_u_prev(pos_counter-1) );
                a_n_d_cur(pos_counter)=( sqrt(1-a*a ) * a_n_u_prev(pos_counter+1)) - ( a * a_n_d_prev(pos_counter+1) );

                %record the probability densities of up,down and all particles
                p_n_u_tmax(p,pos_counter) = (abs(a_n_u_cur(pos_counter)))^2;
                p_n_d_tmax(p,pos_counter) = (abs(a_n_d_cur(pos_counter)))^2; 
                p_n_tmax(p,pos_counter) = (abs(a_n_u_cur(pos_counter)))^2 + (abs(a_n_d_cur(pos_counter)))^2; 
            end;
    %         if(iter == n_iterations)
    %             sigma(t) = ( sum((pos-nmid_rel).^2 .* p_n_tmax) - (sum((pos-nmid_rel).*p_n_tmax)).^2 );
    %             sigma_up(t) = ( sum((pos-nmid_rel).^2 .* p_n_u_tmax) - (sum((pos-nmid_rel).*p_n_u_tmax)).^2 );
    %             sigma_down(t) = ( sum((pos-nmid_rel).^2 .* p_n_d_tmax) - (sum((pos-nmid_rel).*p_n_d_tmax)).^2 );
    %             
    %             avg_pos(t) = sum((pos-nmid_rel).* p_n_tmax);
    %             avg_pos_up(t) = sum((pos-nmid_rel).* p_n_u_tmax);
    %             avg_pos_down(t) = sum((pos-nmid_rel).* p_n_d_tmax);
    %             sigma_count = sigma_count + 1;
    %         end;

            %calculate the sigma 
            p_n(p) = p_n(p) + p_n_tmax(p);
            p_n_u(p) = p_n_u(p) + p_n_u_tmax(p);
            p_n_d(p) = p_n_d(p) + p_n_d_tmax(p);

    %         sig = sig + sigma;
    %         sig_up = sig_up + sigma_up;
    %         sig_down = sig_down + sigma_down;

            %set the a arrays of prev time with values of that of the current
            %time and a arrays of current time with zero values
    %         if(t <tmax)
                a_n_u_prev = a_n_u_cur;
                a_n_d_prev = a_n_d_cur;
                a_n_u_cur = zeros(npos);
                a_n_d_cur = zeros(npos);
    %         if(iter == n_iterations-1) 
    %             p_n_tmax = zeros(1,npos);
    %             p_n_u_tmax = zeros(1,npos);
    %             p_n_d_tmax = zeros(1,npos);
    %         end;

        end;

    end;
end;


p_n = p_n/(n_iterations*tmax);
p_n_u = p_n_u/(n_iterations*tmax);
p_n_d = p_n_d/(n_iterations*tmax);

% sig = sig/(n_iterations*tmax);
% sig_up = sig_up/(n_iterations*tmax);
% sig_down = sig_down/(n_iterations*tmax);

% figure(1)
% plot(pos-nmid_rel,p_n_tmax);
% grid on;

% figure(2)
% plot(T,sigma,T,sigma_up,T,sigma_down);
% grid on;

figure(3)
plot(pos-nmid_rel,p_n(1));
grid on;

% figure(4)
% plot(pos-nmid_rel,p_n_u_tmax);
% grid on;

% figure(5)
% plot(pos-nmid_rel,p_n_u);
% grid on;


% figure(8)
% plot(T,avg_pos);
% grid on;

% figure(9)
% plot(T,sig,T,sig_up,T,sig_down);
% grid on;

